/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"

#include "uhm/matrix/uhm/matrix.hxx"
#include "uhm/matrix/uhm/fla.hxx"

namespace uhm {
  static int qr_flat( int fs, int ss,
                      linal::Flat_ ATL, linal::Flat_ ATR,
                      linal::Flat_ ABL, linal::Flat_ ABR,
                      linal::Flat_ T );
  
  static int qr_hier( int fs, int ss, 
                      linal::Hier_ ATL, linal::Hier_ ATR,
                      linal::Hier_ ABL, linal::Hier_ ABR,
                      linal::Hier_ T );

  void Matrix_FLA_::_qr_create_T() {
    int bmn  = get_hier_blocksize();

#ifdef UHM_HIER_MATRIX_ENABLE

    // ** QR_INC only works for hier matrix
#ifdef UHM_QR_INC_ENABLE
    int t_fs = ( this->fs/bmn > 0 ? ( ( this->fs/bmn + (this->fs%bmn>0) )*bmn ) : this->fs );
    int t_ss = ( this->ss/bmn > 0 ? ( ( this->ss/bmn + (this->ss%bmn>0) )*bmn ) : this->ss );

    flat.T.create_without_buffer(this->datatype, t_fs, t_fs);
    hier.T.create(flat.T, bmn, bmn);
#else
    flat.T.create_without_buffer(this->datatype, this->fs, this->fs);
    hier.T.create(flat.T, bmn, bmn);
#endif


#else
    flat.T.create_without_buffer(this->datatype, this->fs, this->fs);
#endif

    this->create_buffer(UHM_T);
  }
  
  void Matrix_FLA_::qr() {

    // ** create matrix without buffer
    if (!is_created(UHM_T)) 
      this->_qr_create_T();

    // ** decomposition
#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    qr_hier( this->fs, this->ss, 
             this->hier.ATL, this->hier.ATR,
             this->hier.ABL, this->hier.ABR,
             this->hier.T );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix - Level Matrix
    // ----------------------------------------------------------
    qr_flat( this->fs, this->ss, 
             this->flat.ATL, this->flat.ATR,
             this->flat.ABL, this->flat.ABR,
             this->flat.T );
#endif
  }
  
  static inline int qr_flat( int fs, int ss, 
                             linal::Flat_ ATL, linal::Flat_ ATR,
                             linal::Flat_ ABL, linal::Flat_ ABR,
                             linal::Flat_ T ) {
    
    if (fs) 
      FLA_QR_UT( ~ATL, ~T );
    
    if (fs && ss) {

#pragma omp task firstprivate(ATL, ATR, T)
      {
        linal::Flat_ W;

        W.create    ( ATL.get_datatype(), T.get_m(), ATR.get_n() );

        FLA_Apply_Q_UT( FLA_LEFT, FLA_CONJ_TRANSPOSE,
                        FLA_FORWARD, FLA_COLUMNWISE,
                        ~ATL, ~T, ~W, ~ATR ); 
        W.free();
      }

#pragma omp task firstprivate(ATL, ABL)
      FLA_Trsm( FLA_RIGHT, FLA_UPPER_TRIANGULAR,
                FLA_NO_TRANSPOSE, FLA_NONUNIT_DIAG,
                FLA_ONE, ~ATL, ~ABL );

#pragma omp taskwait

      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                FLA_MINUS_ONE, ~ABL, ~ATR, FLA_ONE, ~ABR );
    }

    return true;
  }

  static inline int qr_hier( int fs, int ss, 
                             linal::Hier_ ATL, linal::Hier_ ATR,
                             linal::Hier_ ABL, linal::Hier_ ABR,
                             linal::Hier_ T ) {
    if (fs) {
#ifdef UHM_QR_INC_ENABLE
      linal::qr_inc( ATL, T );
#else
      linal::qr( ATL, T );
#endif
    }

    if (fs && ss) {

#ifdef UHM_QR_INC_ENABLE

#pragma omp task firstprivate(ATL, ATR, T)
      {
	linal::Flat_ W;
	linal::Hier_ W_1;
	
	int bmn   = get_hier_blocksize();
        int t_fs = ( fs/bmn > 0 ? ( ( fs/bmn + (fs%bmn>0) )*bmn ) : fs );

	W.create    ( ATL.get_datatype(), fs, ATR.flat().get_n() );
	W_1.create  ( W, bmn, bmn );
	
	
	linal::apply_q_inc( FLA_LEFT, FLA_CONJ_TRANSPOSE,
				   FLA_FORWARD, FLA_COLUMNWISE,
				   ATL, T, W_1, ATR ); 
	
        W_1.free();
	W.free();
      }

#else

#pragma omp task firstprivate(ATL, ATR, T)
      {
	linal::Flat_ W;
	linal::Hier_ W_1;
	
	int bmn  = get_hier_blocksize();
	
	W.create    ( ATL.get_datatype(), T.flat().get_m(), ATR.flat().get_n() );
	W_1.create  ( W, bmn, bmn );
	
	linal::apply_q( FLA_LEFT, FLA_CONJ_TRANSPOSE,
                               FLA_FORWARD, FLA_COLUMNWISE,
                               ATL, T, W_1, ATR ); 
	
	W_1.free();
	W.free();
      }

#endif

#pragma omp task firstprivate(ATL, ABL)
      {
	linal::trsm( FLA_RIGHT, FLA_UPPER_TRIANGULAR,
			    FLA_NO_TRANSPOSE, FLA_NONUNIT_DIAG,
			    FLA_ONE, ATL, ABL );
      }
      
#pragma omp taskwait
      
      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
			  FLA_MINUS_ONE, ABL, ATR, FLA_ONE, ABR );
    }

#pragma omp taskwait 
    
    return true;
  }
}
